      subroutine grids(qmax, rmax, jump, gridx, nmaxx, nx, 
     >   gridr, nmaxr, nr, jdouble, nmaxj, njdouble)
C
      dimension gridx(nmaxx,3), gridr(nmaxr,3), jdouble(nmaxj)
C
C  This subroutine sets up two grids: gridx and gridr. Gridx will be
C  used to solve the radial Schrodinger equation. Gridr, which is
C  typically every second point of gridx, is used to perform integrations
C  from zero to infinity using Simpson's rule. The structure is such
C  that a wave with momentum QMAX will have NPWAVE points per oscillation.
C  determines HMAX, the largest step. At some point, XDBLE, the intervals,
C  dX, are progressively halved. This is done to accomodate irregular
C  solutions.
C  INPUT:
C   qmax  - the biggest momentum (a.u.) that can be calculated by this mesh.
C   rmax  - the largest "r=x" in the meshes.
C   jump  - gridr(i) = gridx(i*jump)
C   nmaxx - array declaration parameter of first declaration of GRIDX
C   nmaxr - array declaration parameter of first declaration of GRIDR
C   nmaxj - array declaration parameter of first declaration of JDOUBLE
C  OUTPUT:
C   gridx - X grid
C   nx    - Number of X points
C   gridr - R grid
C   nr    - number of R points
C   jdouble - j points where dx doubles
C   njdouble - Number of doublings + 2
C
C  jdouble stores the J index for GRIDX where DX doubles, with the first
C  point set to 1 and the last to NX. Used in the numerov integrations.

C  Set up GRIDX
C  The number of points per oscillation (half-wavelength): NPWAVE
      npwave = 6
C  The number of doubling of intervals is NDOUBLE
      ndouble = 5
      njdouble = ndouble + 2
      if (njdouble.gt.nmaxj) then
         print*,'Increase NMAXJ in call to GRIDS to at least:',njdouble
         stop 'Increase NMAXJ in call to GRIDS'
      end if 
C
C  let NPDBL be the number of points with the same dx per interval
C
      npdbl = 40
C
C  Make sure NPDBL/JUMP is even
C
      npdbl=((npdbl/jump)/2) * 2 * jump
      if (npdbl.lt.4*jump) then
         print*,'Need to have at least 4 equally spaced points in GRIDR'
         stop 'Need to have at least 4 equally spaced points in GRIDR'
      end if 
      hmax = 3.14/float(npwave)/qmax
C
C  2**ndouble * hmin = hmax, therefore
C
      hmin = hmax/float(2**ndouble)
C
C  The value of the X point from which dX is constant, = HMAX, is XDBLE
C  XDBLE = NPDBL * hmin * (2**NDOUBLE-1)
C
      xdble = float(npdbl) * hmin * (2**ndouble-1)
      dx = hmin
      jdouble(1) = 1
      do nd = 2, ndouble + 1
         jdouble(nd) = npdbl * (nd - 1)
      end do 
C
      x = 0.0
      j = 0
      do nd = 1, ndouble
         do nj = 1, npdbl
            j = j + 1
            gridx(j,1) = x + float(nj) * dx
            gridx(j,2) = dx
C
C           Simpson's rule weights
C
            gridx(j,3) = float(mod(j,2) * 2 + 2)
         end do
         x = gridx(j,1)
         gridx(j,3) = 3.0
         dx = dx * 2.0
      end do
      if (abs((dx - hmax)/hmax).gt.1e-4) then
         print*,'DX and HMAX should be equal'
         stop 'DX and HMAX should be equal'
      end if 
      nj = nint((rmax-x)/hmax)
      if (nj+j.gt.nmaxx) then
         nj = nmaxx - j
         print*,'Warning: NX had to be reduced to NMAXX'
      end if 
      do jj = 1, nj
         j = j + 1
         gridx(j,1) = x + float(jj) * dx
         gridx(j,2) = dx
C
C        Simpson's rule weights
C
         gridx(j,3) = float(mod(j,2) * 2 + 2)
      end do
      nx = j
      jdouble(ndouble+2) = j
      print*,'Last X and NX:',gridx(nx,1),nx
C
      nr = nx / jump
      if (nr.gt.nmaxr) then
         nr = nmaxr
         print*, 'Warning: NR had to be reduced to NMAXR'
      end if 
      do i = 1, nr
         j = i * jump
         gridr(i,1) = gridx(j,1)
         gridr(i,2) = float(jump) * gridx(j,2)
         gridr(i,3) = float(mod(i,2) * 2 + 2) * gridr(i,2) / 3.0
      end do 
      do nd = 2, ndouble + 1
         i = jdouble(nd) / jump
         gridr(i,3) = gridr(i,2)
      end do 
      print*,'Last R and NR:',gridr(nr,1),nr
      return
      end
